Spontaneous synchronization of coupled oscillator systems with frequency adaptation 



O 

(N 

D 
tin 



B 



o 

> 

O 
0^ 
O 



X 



Dane Taylor, ' El Edward Ott,^ and Juan G. Restxepo^ 

^Department of Electrical, Computer, and Energy Engineering, 
University of Colorado, Boulder, Colorado 80309, USA 
'^Institute for Research in Electronics and Applied Physics, 
University of Maryland, College Park, Maryland 20742, USA 
^Applied Mathematics Department, University of Colorado, Boulder, Colorado 80309, USA 

(Dated: February 1, 2010) 

We study the synchronization of Kuramoto oscillators with all-to-all coupling in the presence of slow, noisy 
frequency adaptation. In this paper we develop a new model for oscillators which adapt both their phases and 
frequencies. It is found that this model naturally reproduces some observed phenomena that are not qualitatively 
produced by the standard Kuramoto model, such as long waiting times before the synchronization of clapping 
audiences. By assuming a self-consistent steady state solution, we find three stability regimes for the coupling 
constant k, separated by critical points fci and ^2: (i) for k < k\ only the stable incoherent state exists; (ii) for 
k > k2, the incoherent state becomes unstable and only the synchronized state exists; (iii) for k\ < k < k2 
both the incoherent and synchronized states are stable. In the bistable regime spontaneous transitions between 
the incoherent and synchronized states are observed for finite ensembles. These transitions are well described 
as a stochastic process on the order parameter r undergoing fluctuations due to the system's finite size, leading 
to the following conclusions: (a) in the bistable regime, the average waiting time of an incoherent— ^-coherent 
transition can be predicted by using Kramer's escape time formula and grows exponentially with the number of 
oscillators; (b) when the incoherent state is unstable (k > k2), the average waiting time grows logarithmically 
with the number of oscillators. 



I. INTRODUCTION 

Many natural and engineered systems can be described as 
an ensemble of heterogeneous limit-cycle oscillators influenc- 



ing each other Examples include glycolytic oscillations in 
yeast cell populations jTO, pedestrians walking over a bridge 
arrays of Josephson junctions isj], the power grid 101, 



y;east cell populations pedestrians walking over a bridge 

rrays ^ ^ ^ 

lasers and some species of fireflies |6l]. A central issue 
is that of understanding the mechanism of coherent behavior 
that is often observed for these systems. 

The Kuramoto model ^ (for a review, see jstl) addresses 
this problem by considering the simplified case in which os- 
cillators are all-to-all coupled, and each oscillator, labeled by 
an index n, has an intrinsic frequency w„ and an oscillator 
state that can be specified solely by its phase angle 6'„. The 
evolution of the phase of each oscillator n is given by 



k ^ 

1 



(1) 



where N is the number of oscillators, m = 1, 2, . . . , indi- 
cates different oscillators, and fc is a parameter that represents 
the strength of the coupling between oscillators. Kuramoto 
found that for the N ^ 00 coupling limit (approximating the 
typical case of large N arising in most applications), the col- 
lective behavior of the oscillator ensemble, quantified by the 
order parameter r = | X]m=i 6xp(i6',„)|, undergoes a transi- 
tion from incoherence (r = 0) to synchronization (r ~ 1) as 
the coupling strength is increased past a critical value kc- The 
Kuramoto model provides a simple mathematical model cap- 
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turing the essential mechanisms for synchronization of limit- 
cycle oscillators. Despite its long-standing status as a classical 
model of synchronization, some advances in the theoretical 
understanding of Kuramoto-type models have been achieved 
only very recently (e.g., Refs. 121). 

Due to the ubiquity of synchronization phenomena in com- 
plex systems, there is current interest in understanding the ef- 
fect of network structure interactions and adaptation on the 
synchronization of oscillators ifioll . The new model presented 
in this Article aims to investigate synchronization in coupled 
oscillator systems where each oscillator's natural frequency 
U!n slowly adapts, while being subjected to random noise-like 
fluctuations. One motivation for considering adaptive syn- 
chronization is the observation that, in many biological sit- 
uations, synchronization seems to serve a useful function. A 
fairly clear example is the synchronization of pacemaker cells 
in the heart. Another example is the observed evolving pat- 
terns of neuronal synchrony in the brain, which have been 
conjectured to play a key role in organizing brain function. In 
some cases, adaption of frequencies has been experimentally 
observed: fireflies of the species pteroptyx-malaccae slowly 
adapt their flashing frequency in response to the flashes they 
observe Assuming the utility of synchronization in such 
biological cases, it is reasonable that there might be evolu- 
tionary pressure for the development of adaptive mechanisms 
that promote synchronization or maintain it in the presence 
of disruptive influences (e.g., noise). In addition, one could 
imagine technological and social situations where adaptation 
to promote synchronization might be relevant. A familiar so- 
cial example is that of an audience clapping their hands and 
seeking to synchronize ifTTl [T2I1 . Therefore, it is important to 
consider the possibility that various mechanisms might oper- 
ate independently to promote synchronization in a noisy envi- 
ronment. In this paper we introduce and analyze a model of 
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all-to-all coupled phase oscillators with noisy frequency adap- 
tation, where, as seems reasonable in the above-cited biolog- 
ical examples, the coupling of the oscillators' phases occurs 
on a faster timescale than the frequency adaptation dynamics. 

The paper is organized as follows. In Sec.|II]we present and 
analyze our model. Analytical results are derived and numer- 
ically tested in Sees. Ill Al and lll Bl respectively. In Sec.HIIlwe 
study the statistics of spontaneous transitions to the synchro- 
nized state due to finite size effects. In Sec.|lV]we discuss our 
results and their relation with previous work. In Sec. [V] we 
present our conclusions. 



We note that we could have added a noise term to the 9 
equation in However, the effect of such a noise term has 
been akeady studied in the Kuramoto model, and it has been 
found that it shifts the transition to synchronization to larger 
values of the coupling strength k, maintaining the same qual- 
itative behavior. Therefore, for analytical simplicity, we will 
not consider this term. As we will see in our numerical sim- 
ulations, a role analogous to fluctuations in will be played 
by the fluctuations resulting from having a finite number of 
oscillators. 



II. FREQUENCY ADAPTION MODEL 

We consider the classical Kuramoto model supplemented 
with a dynamical equation for the evolution of the oscillators' 
natural frequencies cj„. 



k ^ 

9„ = tJ„ + — ^ sin(6'm - 9n) 

m— 1 

N 



(2) 
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where t is assumed to be much larger than the spread in the 
oscillator period and ?]„ is a Gaussian uncorrelated noise term 
such that {i]n{t)'rim{t')) ~ 2D5nmS{t~t'), where (•) is an en- 
semble average and 5nm is the Kronecker delta. The motiva- 
tion for the different terms in this natural frequency adaption 
term is the following: 

• We assume that each oscillator n may only have knowl- 
edge of the aggregate input to it from the other oscilla- 



tors, fcEm=iSin 



On)- This frequency-coupling 



term was originally introduced in Ref. [Q], who consid- 
ered frequency adaptation without phase coupling. 

• The form of the coupling guarantees that if the phase of 
oscillator n is behind (ahead of) the average phase [so 
that sin(6'„i — dn) is, on average, positive (negative)], 
its frequency increases (decreases). 

• Frequency adaptation occurs on a time scale r, much 
slower than the phase dynamics. 

• The intrinsic frequencies a;„ are subject to random 
noise ?]„. This is partly motivated by observations of 
frequency drift in biological oscillators ifisll . 

• The confining potential Vn{(j~>) represents physical 
mechanisms that, depending on the application, con- 
strain the natural frequencies to some reasonable range. 



2]i is related to that 
1511: however, the 



The dynamics we find for our system 
for the Kuramoto model with inertia il4 
differences are significant and will be discussed in Sec. |IV] 
Also in Sec. |IV] we discuss the relation between (|2]i and a 
model for circadian rhythms that implements wandering, un- 
coupled frequencies Ill6ll . 



A. Model Analysis 

We consider A'^ to be very large and adopt a continuum de- 
scription. Thus, we assume that the ensemble of oscillator 
intrinsic frequencies can be regarded as being drawn from a 
continuous distribution function G(aj, t). 

We analyze our proposed model (|2]i by using the assumed 
separation of timescales between the oscillator dynamics and 
the frequency adaptation. Rewritting Eqs. (|2|l in terms of the 
mean field re*''' = Ylim=i e*^"' we obtain 



LOn - kr sm{dn - V'), 



r ""^fcrsinf 



(3) 
(4) 



Here we have dropped the subscript n on the potential Vn, as 
we will henceforth consider all oscillators to have the same 
Vn- In addition, we will assume that Viyj) = V{~uj). 

Since the frequencies vary on a timescale much longer than 
the phases, on the fast time scale we can approximate w„ in 
Eq. ^ as constant. As we shall soon see, it is relevant to 
assume that G{uj, t) is symmetric in lj and monotonically de- 
creasing away from its maximum. Furthermore, without loss 
of generality, it suffices to take the maximum of G to occur at 
w = (if it does not, it can be shifted to by the change of 
variables uj uj+il, 9 9 — fit). As originally noted by Ku- 
ramoto, in the saturated state (i.e., r constant on the fast time 
scale), the phase dynamics is of two types depending on the 
value of ujn- For |a;„| < kr oscillator n is said to be "locked" 
and its phase settles at a value given by 



sm 



-0) = iOn I {kr). 



(5) 



For |aj„| > kr the phase is said to "drift" and 6'„ continually 
increases (decreases) with time for cj„ > kr {ujn < ^kr). For 
a given frequency |a;| > kr, the drifting oscillators have a dis- 
tribution of phases p{9,uj) determined from the conservation 
of oscillator density by the condition pd9/dt = constant. This 
yields 



p{9,uj) 



[krf 



2-k\io — krsvvJyf) ~ -0)1 ' 



(6) 



where the factor lo'^ — [kr)'^ /{2'k) normaUzes p{9,uj) so 
that /^^ pd9 = l. 

Still invoking the time scale separation, and consequently 
assuming that the deterministic term in Eq. ^ can be aver- 
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aged over time, we obtain an approximation to Eq. (|4| for the 
drifting oscillators 



p{9, uj„) sm{9 — i]j)d9 + rjn — dV/dujn- 



(7) 

[Here the deterministic term has been replaced by its time av- 
erage, while the fast-varying noise term has been retained. 
This can be justified by noting that the difference between the 
original equation, Eq. dU, and the equation where the deter- 
ministic part was averaged, Eq. (|7]i, is what would be obtained 
in the noiseless case. Since in the noiseless case this differ- 
ence can be argued to be small by averaging, we conclude our 
procedure is justified.] 

Integrating Eq. (|7J and recalling that for entrained oscilla- 
tors kr sin(6'„ — V-') = ^n, Eq. ^ can be rewritten as 



where 



UJn = h{uJn) +r]n- dV/diOn, 



(8) 



< kr. 



-uj/t + sign(a;)-\/aj2 — (kr)"^ /t, > kr. 



We now seek a steady state solution (in a statistical sense) 
for Eqs. (|2]i. More precisely, for a given value of fc, we 
seek a time-independent probability distribution of frequen- 
cies G and a value of r that make Eqs. (|2]i consistent. Such a 
steady-state frequency distribution can be obtained by solving 
the time-independent Fokker-Planck equation corresponding 
to Eq. (Ull 



d 

duj 



duj 



= D- 



doj^ ■ 



(9) 



The solution of Eq. (|9]l with no-flux boundary conditions 

(dG/dtu = at either uj ~ ±oo or w = ±L) is G cx 
exp[/ h{uj)du} / D—V{uj) / D], from which we obtain Eq. ( fTOb . 
where ct^ = Dt: 



G{uj, kr) cx 



M(l-Vl-(fcr/'^) 



exp 



\u\ < kr, 

^{l- v/1 - {krH^) - V{u)/a^] , \oj\ > kr. 



(10) 



In all our numerical plots we use D = 0.01 and t = 50, 
which yields cr^ = 1/2. 

This distribution depends on the value of the order parame- 
ter r. In order to make this solution self-consistent, the value 
of r has to be determined from the classical Kuramoto results 
corresponding to an ensemble of oscillators with frequency 
distribution G{lu, kr). That is, r is equal to the average of 
exp(i6') over all the oscillators. As shown, e.g., in Refs. lU, 
the average over drifting oscillators (|i^| > kr) is zero, and r is 
thus determined entirely by the locked oscillators (|a;| < kr) 
whose phase angles are given by Eq. Thus we obtain 

/kr 
G{uj, kr)^/l - {uj/kr)^duj. 
-kr 

Besides the solution r = 0, other possible values of r are 
given by the solutions of the nonlinear equation [jst] 



l = fcy G{zkr, fcr)Vl - z^dz. (11) 

We now study numerically the solutions of Eqs. ( fTOb and 
(fTTT i. As we will later argue, noise in the 6 evolution (either 
extrinsic or due to finite N) causes the dynamics to be insensi- 
tive to the choice of confining potential V{uj). Therefore, for 
simplicity, we will choose an infinite potential well to simplify 



the analysis. The potential is defined as V{Ld) = if |a;| < L, 
and V{uj) = oo otherwise. This corresponds to frequencies 
that evolve freely in a box of size 2L with hard walls. We 
will later show that L can be chosen large enough so that the 
dynamics is insensitive to its value. Except for Fig.[Tl all our 
plots use L = 5. 

Solving Eq. ( fTTT i numerically, a bifurcating pair of solutions 
Ts (fc) and ru{k) is found to appear at a finite value of the cou- 
pling strength A: = fci as shown in Fig. [T] The upper branch 
(solid black line) in these figures was numerically found to 
be stable, while the lower branch (colored lines in Fig. [T] and 
dashed line in Figs. [SJi andO was numerically found to be 
unstable. The trivial solution r was numerically found to 
be stable until the lower branch crosses ?- = at a value of 
the coupling constant k = k2 (see Fig.[T]). For larger coupling 
strength, k > k2, the nonzero unstable solution disappears 
and the solution ?- = becomes unstable. 

Thus, three regimes are found with our model. For k < ki, 
the oscillator ensemble is incoherent, as only the r = stable 
solution exists. For k > k2, only the synchronized state is sta- 
ble. These correspond to the traditional regimes of the origi- 
nal Kuramoto model without adaptation Jsll • The third regime 
corresponds to intermediate coupling strengths on the finite 
interval fci < k < k2, where both the synchronized and inco- 
herent states are locally stable, and whose basins of attraction 
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FIG. 1. Stable (upper black line, rs(fc)) and unstable (lower colored 
lines, ru{k)) branches are shown for D = 0.01, r — 50, and L — 
5, 10, 15, 20 and the curve kr — y/2a (dotted line), above which the 
frequency distribution is normalizable. The values of ki and k2 for 
L = 5 are indicated by vertical arrows at ki « 1.8 and k2 « 6.37. 



are separated by the unstable solution of Eq. (fTTl i. A similar 
regime was also found in Ref. lITill for an inertial version of 
the Kuramoto model without noise [in which 9n is replaced 
by m9n + On in O]- In Sec. |III]we address the important 
issue of noise-induced spontaneous transitions between stable 
solutions which was not addressed in Ref. ifTill . 

To study the behavior of our model close to the incoherent 
state r = 0, we expand G in Eq. ( fTOl ) for {kr/ut)'^ <^ 1, find- 
ing that the oj dependence of G for large \uj\ is G{uj, kr) ^ 
{kr/uiY^^^ 1'^'^"^ . Thus the frequency distribution G(aj) can 
be normalized in (— oo, oo) if kr > \/2a. We conclude that, 
as long as kr > V2a, G should be insensitive to L if the bulk 
of G is contained within {—L,L). In contrast, for kr < \f2G, 
the frequency distribution is not normalizable in (— oo, oo) 
and thus we expect G to be broadly distributed in (— L,i), 
and the dynamics to depend on the value of L. In particular, 
the distribution of frequencies for the incoherent state r = is 
uniform with G(uj) ~ l/(2i). More generally, G should be 
insensitive to the specific form of the confining potential V as 
long as V is negligible in the synchronized state. Via) <C cP' . 
This can be interpreted as requiring that V does not itself pro- 
mote synchronization. 

Figure [T] shows the stable and unstable branches for var- 
ious values of L and the curve kr = \/2g (dotted line), 
above which the frequency distribution is normalizable. Note 
that the stable solution rs(uj) to Eq. ( fTTT i is above the curve 
kr = \f2a in Fig. [T] and is, as expected, insensitive on the 
value of L. One interesting result from Fig. [T]is that the upper 
critical coupling strength k^ depends on the frequency bound 
L. Recalling that G(y:) = 1/ (2i) for the incoherent state, one 
can integrate Eq. ( fTTT i to find fc2 = 4i/ tt. In reality, physical 
limitations typically bound the natural frequency distribution. 
However, it is also interesting to consider the unbound case 
(y = 0, or L = oo) which leads to the following scenario: the 
oscillators' natural frequencies wander in (— oo, oo); fc2 = oo 
and the incoherent state r = remains stable for all coupling 
strengths. However, it is important to note that, even though 
the incoherent state for L — oo remains stable for arbitrarily 
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FIG. 2. (a) Order parameter r obtained from numerical simulation 
of Eqs. l[2j for decreasing values of k for A'^ = 10* (triangles) with 
D — 0.01, T = 50, and L — 5. The solid and dashed lines indicate 
stable rs{k) and unstable ru{k) solutions to Eq. i ll It . respectively. 
The letters a, b and c indicate values of k at which the frequency 
distribution is sampled for Fig. [2{). (6) Frequency distribution ob- 
tained directly from Eqs. ^ (symbols) and from Eq. l llOt (black 
lines). Curves labeled a, b, and c correspond to = 1, 2, and 4, 
indicated by arrows in Fig. |2^. The dashed vertical lines indicate 
±kr for k — 2. 



large k, its basin of attraction shifts to zero as fc oo (the 
unstable equilibrium r„ approaches zero as /c — >■ oo). The 
situation is analogous to that applying in the study of the tran- 
sition of stable laminar pipe flow to turbulence (e.g., Ref. ITtIi 
and references therein). As in pipe flow, this situation points 
to the possibly crucial role of noise which we address subse- 
quently. 



B. Model Simulation 

In order to test our theoretical results, we compared them 
with direct simulation of Eqs. Q using t = 50 and D = 0.01. 
Due to the stability characteristics of the solutions, hystere- 
sis phenomena and dependence on the initial conditions are 
expected. To probe these characteristics, we let L = 5 and 
initiate a simulation with strong coupling (fc > k2) and with 
the phases and frequencies of the oscillators clustered around 
9n ~ and cj„ w 0. The oscillators remain synchronized, 
and their natural frequencies adopt a distribution given by 
Eq. ( [Tol l. For a given value of k, we simulate Eqs. (|2| for 
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1000 seconds and then decrease the value of fc by 0.1, keep- 
ing the values of the phases and frequencies (this corresponds 
to a coarse grained rate dk/dt w 10^^). As this process is 
repeated and the value of k decreases below ki, the synchro- 
nized solution disappears and the oscillators desynchronize. 
FigurellJi shows the value of r obtained by this process (trian- 
gles). The solid and dashed lines indicate the stable rs(fc) and 
unstable ru(fc) solutions, respectively obtained from Eq. ( fTTT i. 
The numerically obtained values of r follow the stable branch 
found theoretically. 

In Fig. |2f) we show the steady-state frequency distribution 
observed at values of k corresponding to the arrows labeled a, 
b, and c in Fig.|2}i. The black solid, dashed, and dotted lines 
indicate the theoretical expression given by Eq. ( fTOl i normal- 
ized on w e [—5, 5] for cases a, b, and c. The cross, circle, and 
square symbols show the corresponding observed frequency 
distributions which are in good agreement with the theory. 

To observe hysteresis phenomena similar to that noted 
in lfl4ll . the system was was brought to steady state with 
a dispersed frequency distribution described by Eq. (fTOl i 
for small coupling strength (fc < fci). The coupling 
strength k was slowly increased until the system underwent 
an incoherent ^-synchronized state transition at the transition 
coupling strength k*, which is found on the interval fci < 
k* < k2- The precise value of k* fluctuates slightly from run 
to run, but its mean is observed to depend on the ensemble 
size N. This is shown in Fig. [3] where k* approaches fc2 as 
N increases, as previously noted in lITill for the inertial Ku- 
ramoto model. This is due to fluctuations in the order param- 
eter Ar ^ N^^/^ resulting from the system's finite size: it 
is hypothesized that fluctuations cause the system to cross the 
barrier imposed by the unstable solution to Eq. ( fTTT ) (dashed 
line in Fig. [3]). When the size of these fluctuations becomes 
large enough to place r above the unstable solution, the oscil- 
lators begin to synchronize and the value of the order parame- 
ter increases to the value corresponding to the stable solution 
(upper solid line). 

It should also be noted that, for this simulation with tempo- 
rally increasing coupling strength, the k* approach fci as the 
simulation duration for each fc is increased. In other words, 
the hysteretic nature of this system depends not only on the 
size of the network (as noted in 13), but also on the rate at 
which the coupling strength fc is varied. We hypothesize this 
phenomenology to also describe other Kuramoto-type models 
with hysteretic behavior (e.g., ifTii ITsIl ). The fluctuations of 
the order parameter r are stochastic, and thus the time required 
for the transition to occur is a random variable. The longer a 
simulation is run at constant coupling strength fci < fc < fc2, 
the more likely an incoherent — > synchronized transition has 
occurred. In fact, oscillations between states, as hypothesized 
in ifTill . were observed for our model in this bistable regime 
(see Fig. |7]i. Describing such spontaneous state transitions is 
the focus of the next section of our paper 
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FIG. 3. For increasing coupling strength, synchronization occurs for 
each network when the order parameter fluctuations Ar allow r to 
surmount the barrier of the unstable solution r„(fc) (dashed line). 
Simulation used D = 0.01, r = 50, and L = 5. Note that the transi- 
tion coupling strengths k* approach k2 as network size A'^ increases. 

III. SPONTANEOUS STATE TRANSITIONS 

Given the observed phenomenology of fluctuations driving 
the system from one stable solution to another across an un- 
stable solution, it is natural to conjecture that, for a fixed value 
of fc, the average time Tsync {N) for a transition from the inco- 
herent state r « to the coherent state r ^ 1 can be obtained 
by treating the problem as an escape over a potential barrier 
under the influence of random noise (see Fig. |4]i. Conceptu- 
ally, it is helpful to relate such transitions to a Brownian par- 
ticle moving from one equilibrium to a second by traversing 
an energy barrier under the influence of random noise. For the 
case of oscillator system state transitions, fluctuations of the 
order parameter Ar occur due to a network's finite size and 
are akin to random noise. In addition, in some applications, 
Eq. (|2]i may be subject to extrinsic noise iU. 

For the traditional Kuramoto model, understanding finite 
size fluctuations Ar has been a major area of interest ll8l [l9ll . 
In general, fluctuations are typically 0{N~^/^), although it 
has been shown that these fluctuations increase in amplitude 
near the critical coupling fcc lfl9ll for a traditional Kuramoto 
oscillator system. Similarly, for our model, fluctuations in 
r were observed to be larger in the bistable regime than in 
the traditional Kuramoto regimes. However, as with the tradi- 
tional Kuramoto model, further study of these fluctuations for 
our model remains open to future research. 

A. State Transition Analysis 

In order to study the statistics of spontaneous synchroniza- 
tion transitions, we will assume that finite-size fluctuations 
can be described approximately as produced by uncorrected 
Gaussian noise acting on the 1 -dimensional dynamics of the 
order parameter Treating finite-size fluctuations as an un- 
corrected Gaussian noise term has already proven sucessful 
in studying synchronization of Kuramoto oscillators in net- 
works 12011 . Consequently, let us assume that the macroscopic 
dynamics of the order parameter r can be described by a 



FIG. 4. State transitions parameterized by r for fci < fc < k2 are 
schematically shown as Brownian motion in a 1 -dimensional energy 
landscape with two stable equilibria. 

Langevin equation of the form 

f = -C/'(r,fc)+i(0, (12) 

where U{r,k) is an unknown pseudo-potential, U'{r,k) ~ 
dU jdr, and L(t) is an uncorrelated Gaussian noise term such 
that {L{t)) ^ and {L{t)L{t')) = 2T5{t - t'). Since the 
noise represents finite-size fluctuations, the diffusion coeffi- 
cient r will be assumed to be inversely proportional to N, 
or r oc Note that this is consistent with Ar being 

0{N^'^/'^) for the dynamics of r modeled as a linear Ornstein- 
Uhlenbeck process for the incoherent state with k < ki. 

In the bistable regime, ki < k < k2, we assume U{r, k) to 
be of the form shown in Fig. HI Potentials of this type have re- 
ceived much attention in the literature for studying Brownian 
motion in bistable potentials and for describing chemical reac- 
tions. We will draw on this research and use Kramer's escape 
time equation |l2l!] , which describes the mean first-passage 
time Tesc for a particle subject to random noise with diffu- 
sion coefficient F to escape over a potential barrier of height 
h, and is given by log(resc) oc h/T. Recalling that F a 
we conclude that the mean first-passage time (i. e., wait time 
before synchronization) for our bistable Kuramoto system de- 
pends exponentially on N, yielding Tgync cc e^^ for some 
constant K. 

A similar analysis can also be done on the regime where 
the incoherent state is unstable, where we are interested in 
the average time required for an incoherent system (r ^ 0) 
to synchronize. To first order, the dynamics for small r is 
described by f = ar + L{t), with a being a positive constant. 
Taking r(0) = and setting {r{t*Y) = r*^, we can estimate 
for large N the time t* it takes for the order parameter to reach 
a given threshold r = r* > ^T/a as t* ~ log F^^ ~ log N. 
Thus, the waiting time Tgync grows logarithmically with N in 
the strong coupling regime (k > fc^ in the Kuramoto model or 
fc > fc2 in our model). 

Although this paper focuses on the model described by 
Eqs. (|2]i, the above estimates may apply to other Kuramoto- 
type models |[ll[il[il. 

B. State Transition Simulation 

To test the previous findings, statistics were compiled for 
our adaptive Kuramoto system by simulating 100 realizations 
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FIG. 5. Synchronization time Tsync averaged over 100 realizations 
as a function of the number of oscillators A'^ for fc = 6, which is 
within the bistable regime. (D = 0.01, r = 50, and L = 5) 
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FIG. 6. Synchronization time Tsync averaged over 100 realizations 
as a function of the number of oscillators A*' for k = 7 > k2. (D = 
0.01, T — 50, and L — 5). Note that the scale is different than that 
ofFig.S 

of synchronization for an initially incoherent system. For each 
realization, at a constant coupling strength k the initial natural 
frequencies and phases were chosen randomly (6'„ uniform in 
[0, 27r), and a;„ uniform in [—5,5]). Once the order parameter 
exceeded a given threshold r* ensuring synchronization had 
occurred, the time before synchronization was recorded and 
simulation stopped. 

Statistics of incoherence^synchronization transitions for 
the bistable regime are shown in Fig.|5] where log(Tsj,„c(^)) 
vs. N is plotted for k = 6. Tgync is defined as the average 
time required for the order parameter to first reach r* ~ 0.7. 
In principle any coupling strength k within the bistable regime 
could be used; however, to decrease simulation time k was 
chosen to be close to fc2 ~ 6.37. Error bars are included to 
show statistical uncertainty. As the plot shows, log{Tsync{N)) 
is well described by a straight line, which is consistent with 
the supposition that the transition times can be described by 
Kramer's escape time formula. 

For comparison, Tgijnc is shown in Fig. |6]for synchroniza- 
tion with the incoherent state being unstable {k > k2)- From 
this figure we confirm that Tgync cc log N, which is consis- 
tent with unstable exponential growth of perturbations from 
the r ~ incoherent state. 

Figure |7] shows fluctuations between the synchronized and 
incoherent states for a case where the coupling strength is 
within the bistable range. Note that since transitions between 
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FIG. 7. Spontaneous bidirectional transitions between the synchro- 
nized (dashed) and incoherent (dotted) states are observed for — 
10, k = 1.9, D = 0.01, r = 50, and L = 5. Note that because 
of the small system size, the incoherent state has an average order 
parameter of (r) ~ 0.4. 



States are related to the height h of the pseudo-potential bar- 
rier relative to each respective equilibrium (see Fig. |4|, fluc- 
tuations between states can only be observed when the barrier 
heights are roughly equal and when the system is observed 
for a duration in which transition-events should occur For 
example, if the barrier height is large and the finite system 
is large (large N), the order parameter r will undergo small 
fluctuations and state transitions would be rare. At the same 
time, if the barrier height is much larger for a particular state, 
then the system will remain in that state for the majority of 
time and transitioning out of that state would also be rare. For 
the model parameters chosen in our simulation, we found that 
spontaneous bidirectional transitions could only be observed 
for small numbers of oscillators {N = 10 in Fig. |7]i and for 
coupling strengths in the bistable regime just above ki (be- 
low which the coherent solution disappears). In general, for 
ki < k < k2, we find that synchronized— >incoherent tran- 
sitions are very rare, implying that the barrier height for the 
synchronized state is generally larger than the barrier height 
for the coherent state (as shown schematically in FigSJi. 



IV. DISCUSSION 

Our results discussed above are in striking agreement with 
observations of rhythmically clapping audiences ifTTl [l2ll . In 
particular, as opposed to the behavior of the classical Ku- 
ramoto model without adaption, the transition to synchronized 
clapping occurs after a relatively long waiting time, and once 
it starts the order parameter quickly achieves its steady state. 
Previous models of this phenomenon have artificially altered 
the frequency distribution ifllll or introduced additional dy- 
namics such as a time-dependent tendency of the oscillators 
to synchronize 1 12|. In contrast, the long waiting times arise 
in our model as a natural consequence of the dynamics. Al- 
though we have found that all-to-all coupling leads to waiting 
times that depend exponentially on the number of oscillators, 
shorter waiting times are expected for local coupling such as 
that describing clapping synchronization in a large venue. 

Another possible application of our model is circadian 



rhythms II13I1 . which have been modeled by ensembles of Ku- 
ramoto oscillators with drifting, nonadaptive frequencies lfl6ll . 
Because of the importance of synchronization in this system, 
evolutionary pressures might have led to frequency adapta- 
tion. Our model generalizes previous models jlql by allowing 
for frequency adaptation. By removing frequency coupling 
(i.e., T oo) and assuming a quadratic form for the poten- 
tial V{uj), our model [Eqs. (|2|i] recovers the model of coupled 
circadian oscillators presented in lfl6ll . 

Our results are somewhat related to the Kuramoto model 
with inertia (Eq. ([TJ with 9n replaced by m9n+6n flE 
which is equivalent to 



N 



dVn{u}n)/du}n + ^ siii(6',„ - 6'„) 



(13) 



Vn 



m— 1 



where y„(ci;„) ~ ^(a;„ — 51„)^, with rt„ constant for each os- 
cillator. However, the differences between this model and our 
model are significant. First, in contrast to (fT3T l. our model cou- 
ples both phases and frequencies. Second, as a consequence of 
our two types of coupling, we are able to introduce two time 
scales, with the frequency adaptive time scale being slower 
than that of the phase dynamics. We believe that this two time 
scale dynamics will be crucial to the modeling of the various 
potential applications mentioned in Sec.|I](e.g., clapping audi- 
ences). [Note that simulations were conducted to investigate 
the effect of closing the timescale gap. By keeping cr constant 
and reducing t (i.e., by also increasing D), it was found that 
no qualitative differences were observed as long as t > 5.] 

The analysis presented here to describe fluctuation-induced 
spontaneous transitions from incoherence to synchronization 
for our adaptive model could also be applicable to other 
Kuramoto-type systems with hysteretic behavior. Such sys- 
tems include Kuramoto models with an added inertial term 
ifTii [Tsll and situations where there is a heterogeneous dis- 
tribution of interaction time delays 1 1 8i1 . Various questions 
remain to fully understand the dynamics of the observed 
transitions. While order parameter fluctuations are typically 
0{N~^^'^), this is not always the case and a better understand- 
ing of these fluctuations is needed. Similarly, the existence of 
a pseudo-potential U (r, k) was assumed [Eq.[T2l, but its shape 
and dependence on k remain to be investigated. 



V. CONCLUSIONS 

We have presented a new model to study the synchroniza- 
tion of Kuramoto oscillators that are able to slowly adapt their 
natural frequencies to promote synchronization, but are in- 
hibited from doing so completely by the influence of noise. 
We found that the interplay of noise and adaptation results in 
bistability and hysteresis. In the bistable regime, finite size 
effects induce incoherent^synchronized state transitions (or, 
when N is small, vice versa), which are well described as a 
1 -dimensional Kramer escape process on the order parameter 
r. For an oscillator ensemble governed by our adaptive model 
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with all-to-all coupling, it was shown that the time Tsync re- 
quired for the system's state to transition from incoherent to 
synchronized depended exponentially on N in the bistable 
regime (fci < k < and logarithmically for strong cou- 
pling {k> k2). 

To our knowledge, this work is the first to analyze spon- 
taneous synchronization at constant coupling strength as a 



1-dimensional stochastic escape process. It is expected that 
the analysis presented in this paper is also valid for other 
Kuramoto-type models with hysteretic behavior lfl4l [Tsi [Tsll . 
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